Load the data

library(ncdf4)
nc <- nc_open("wind.change.hgt.ivg.nc")
spd_delta <- ncvar_get(nc, "spd")
spd95_delta <- ncvar_get(nc, "spd95")
Lat <- ncvar_get(nc, "lat")
Lon <- ncvar_get(nc, "lon")

ncFiles <- list.files(pattern = "*.nc")
nc <- lapply(ncFiles[-3], nc_open)

spd_p <- ncvar_get(nc[[3]], "spd")
spd_f <- ncvar_get(nc[[4]], "spd")
spd95_p <- ncvar_get(nc[[1]], "spd95")
spd95_f <- ncvar_get(nc[[2]], "spd95")

Plotting rounties

# Produce map in changes (red increase, blue decrease, white no change)
two.cols <- function(dat, nlevels, col = c("blue", "white", "red")){
  n = nlevels
  z = c(dat)
  n1 = floor(n * (abs(range(z)[1])) / (range(z)[2] - range(z)[1]))
  n2 = n - n1
  col1 = colorRampPalette(col[1:2])(n1)
  col2 = colorRampPalette(col[2:3])(n2)
  cols = c(col1, col2)
  return(cols)
}

plot_WRF_wind <- function(lon = lon, lat = lat, dat, col = "gray",
                          image.col = tim.colors(), title,...){
  par(mar = c(0, 4, 1, 1))
  maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
            mar = c(0, 3, 1, 1))
  #axis(1); axis(2)
  image.plot(lon, lat, dat, add = T, horizontal = T, col = image.col,
             ...)
  maps::map("world", xlim = c(-145, -45), ylim = c(20, 60), add = T,
            col = col)
  maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
            add = T, col = col)
  title(title)
}

Selected locations

select <- matrix(c(-120.6, 33.8017, -119.735, 33.4773,
                   -116.02, 43.2696, -114.69, 43.7076,
                   -124.419, 44.2998, -97.5567, 32.0997, 
                   -99.8971, 48.2343, -82.3954, 36.1412,
                   -73.0984, 40.606, -82.4075, 41.6163
                   ), 10, 2, byrow = T)

## Restrict the spatial domain
xid <- 125:550; yid <- 30:330
lon <- Lon[xid, yid]; lat <- Lat[xid, yid]
## Search gird cell indices for the selected grids
library(fields)
dist2select <- rdist.earth(select, cbind(c(lon), c(lat)), miles = F)
id <- apply(dist2select, 1, which.min)
index <- function(x, row = 426) matrix(c(x %% row, x %/% row + 1), 1, 2)
ID <- sapply(id, index)

Annual Mean Wind Speed

# Present
plot_WRF_wind(lon, lat, dat = spd_p[xid, yid], title = "1995-2004 mean wind speed (m/s)")
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Future
plot_WRF_wind(lon, lat, dat = spd_f[xid, yid], title = "2085-2094 mean wind speed (m/s)")
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Change
delta <- spd_f[xid, yid] - spd_p[xid, yid]
plot_WRF_wind(lon, lat, dat = delta,
              title = "Change in mean wind speed (m/s)",
              image.col = two.cols(delta, 200))
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Ratio
ratio <- log(spd_f[xid, yid] / spd_p[xid, yid])
ticks <- seq(-0.2, 0.1, 0.05)
plot_WRF_wind(lon, lat, dat = ratio,
              title = "Ratio of mean wind speed",
              image.col = two.cols(ratio, 200),
              axis.args = list(at = ticks,
                               labels = round(exp(ticks) - 1, 2)))
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

95th Percentile Wind Speed

# Present
plot_WRF_wind(lon, lat, dat = spd95_p[xid, yid], title = "1995-2004 95th percentile wind speed (m/s)")
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Future
plot_WRF_wind(lon, lat, dat = spd95_f[xid, yid], title = "2085-2094 95th percentile wind speed (m/s)")
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Change
delta <- spd95_f[xid, yid] - spd95_p[xid, yid]
plot_WRF_wind(lon, lat, dat = delta,
              title = "Change in 95th percentile wind speed (m/s)",
              image.col = two.cols(delta, 200))
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

# Ratio
ratio <- log(spd95_f[xid, yid] / spd95_p[xid, yid])
ticks <- seq(-0.2, 0.1, 0.05)
plot_WRF_wind(lon, lat, dat = ratio,
              title = "Ratio of 95th percentile wind speed",
              image.col = two.cols(ratio, 200),
              axis.args = list(at = ticks,
                               labels = round(exp(ticks) - 1, 2)))
for (i in 1:10) points(lon[ID[1, i], ID[2, i]], lat[ID[1, i], ID[2, i]],
       pch = "+")

Semi-Automated Search for other Locations

## Large increases in 95th precentile

temp <- sapply(order(spd95_delta[xid, yid], decreasing = T)[1:600], index)
maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "red", cex = 0.5)

## Large decreases in 95th precentile

temp <- sapply(order(spd95_delta[xid, yid])[1:600], index)
maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "blue", cex = 0.5)

## Large increases in mean

temp <- sapply(order(spd_delta[xid, yid], decreasing = T)[1:600], index)
maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "red", cex = 0.5)

## Large decreases in mean

temp <- sapply(order(spd_delta[xid, yid])[1:600], index)
maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "blue", cex = 0.5)

## Large increases in 95th precentile/mean ratio
temp <- sapply(order(log(spd95_f[xid, yid] / spd95_p[xid, yid]) - log(spd_f[xid, yid] / spd_p[xid, yid]), decreasing = T)[1:600], index)

maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "red", cex = 0.5)

## Large decreases in 95th precentile/mean ratio
temp <- sapply(order(log(spd95_f[xid, yid] / spd95_p[xid, yid]) - log(spd_f[xid, yid] / spd_p[xid, yid]))[1:600], index)

maps::map("world", xlim = c(-145, -45), ylim = c(20, 60),
          mar = rep(0, 4))
maps::map("state", xlim = c(-145, -45), ylim = c(20, 60),
          add = T)
for (i in 1:600)
points(lon[temp[1, i], temp[2, i]], lat[temp[1, i], temp[2, i]],
       pch = "+", col = "blue", cex = 0.5)